Goal: can machine learning methods help us to associate metabolites with leaf length?

Previously (script 11b2) I filtered out unnamed metabolites. Here I keep them all. but remove blank samples

library(glmnet)
Loading required package: Matrix
Loaded glmnet 4.1
library(relaimpo)
Loading required package: MASS
Loading required package: boot
Loading required package: survey
Loading required package: grid
Loading required package: survival

Attaching package: ‘survival’

The following object is masked from ‘package:boot’:

    aml


Attaching package: ‘survey’

The following object is masked from ‘package:graphics’:

    dotchart

Loading required package: mitools
This is the global version of package relaimpo.

If you are a non-US user, a version with the interesting additional metric pmvd is available

from Ulrike Groempings web site at prof.beuth-hochschule.de/groemping.
library(tidyverse)
Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ───────────────────────────────── tidyverse 1.3.0 ──
✓ ggplot2 3.3.3     ✓ purrr   0.3.4
✓ tibble  3.0.4     ✓ dplyr   1.0.2
✓ tidyr   1.1.2     ✓ stringr 1.4.0
✓ readr   1.4.0     ✓ forcats 0.5.0
── Conflicts ──────────────────────────────────── tidyverse_conflicts() ──
x tidyr::expand() masks Matrix::expand()
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()
x tidyr::pack()   masks Matrix::pack()
x dplyr::select() masks MASS::select()
x tidyr::unpack() masks Matrix::unpack()

get leaflength data

leaflength <- read_csv("../../plant/output/leaf_lengths_metabolite.csv") %>%
  mutate(pot=str_pad(pot, width=3, pad="0"),
         sampleID=str_c("wyo", genotype, pot, sep="_"), 
         group=str_c(soil,genotype,trt,sep="_")) %>%
  select(sampleID, group, leaf_avg_std)  %>%
  filter(!str_detect(group, "BLANK"))

── Column specification ──────────────────────────────────────────────────
cols(
  pot = col_double(),
  soil = col_character(),
  genotype = col_character(),
  trt = col_character(),
  leaf_avg = col_double(),
  leaf_avg_std = col_double()
)
length(unique(leaflength$group))
[1] 4
leaflength %>% arrange(sampleID)

get and wrangle metabolite data

met_raw <-read_csv("../input/metabolites_set1.csv")

── Column specification ──────────────────────────────────────────────────
cols(
  .default = col_double(),
  tissue = col_character(),
  soil = col_character(),
  genotype = col_character(),
  autoclave = col_character(),
  time_point = col_character(),
  concatenate = col_character()
)
ℹ Use `spec()` for the full column specifications.
met <- met_raw %>% 
  mutate(pot=str_pad(pot, width = 3, pad = "0")) %>%
  mutate(sampleID=str_c("wyo", genotype, pot, sep="_")) %>%
  filter(soil!="BLANK") %>%
  select(sampleID, genotype, tissue, sample_mass = `sample_mass mg`, !submission_number:concatenate) %>%
  pivot_longer(!sampleID:sample_mass, names_to = "metabolite", values_to = "met_amount") %>%
  
  #adjust by sample mass
  mutate(met_per_mg=met_amount/sample_mass) %>%
  
  #scale and center
  group_by(metabolite, genotype, tissue) %>%
  mutate(met_per_mg=scale(met_per_mg),
         met_amt=scale(met_amount)
  ) %>% 
  pivot_wider(id_cols = sampleID, 
              names_from = c(tissue, metabolite), 
              values_from = starts_with("met_"),
              names_sep = "_")

met

split this into two data frames, one normalized by tissue amount and one not.

met_per_mg <- met %>% select(sampleID,  starts_with("met_per_mg")) %>%
  as.data.frame() %>% column_to_rownames("sampleID")
met_amt <- met %>% select(sampleID,  starts_with("met_amt")) %>%
  as.data.frame() %>% column_to_rownames("sampleID")

get leaf data order to match

leaflength <- leaflength[match(met$sampleID, leaflength$sampleID),]
leaflength

Calc PCAs:

normalized

leaf

met_per_mg.leaf_PCA <- met_per_mg %>% 
  select(matches("_leaf_")) %>%
  prcomp(center = FALSE, scale. = FALSE) #already centered and scaled
names(met_per_mg.leaf_PCA)
[1] "sdev"     "rotation" "center"   "scale"    "x"       
tibble(variance=met_per_mg.leaf_PCA$sdev^2, PC=str_c("PC", 
                                                      str_pad(1:length(met_per_mg.leaf_PCA$sdev), width = 2, pad="0"))) %>%
  mutate(percent_var=100*variance/sum(variance),  
         cumulative_var=cumsum(percent_var)) %>%
  magrittr::extract(1:15,) %>%
  ggplot(aes(x=PC, y=percent_var)) +
  geom_col(fill="skyblue") + 
  geom_line(aes(y=cumulative_var), group="") +
  ggtitle("percent variance explained, named, normalized leaf metabolites")

root

met_per_mg.root_PCA <- met_per_mg %>% 
  select(matches("_root_")) %>%
  prcomp(center = FALSE, scale. = FALSE) #already centered and scaled
names(met_per_mg.root_PCA)
[1] "sdev"     "rotation" "center"   "scale"    "x"       
tibble(variance=met_per_mg.root_PCA$sdev^2, PC=str_c("PC", 
                                                      str_pad(1:length(met_per_mg.root_PCA$sdev), width = 2, pad="0"))) %>%
  mutate(percent_var=100*variance/sum(variance),  
         cumulative_var=cumsum(percent_var)) %>%
  magrittr::extract(1:15,) %>%
  ggplot(aes(x=PC, y=percent_var)) +
  geom_col(fill="skyblue") + 
  geom_line(aes(y=cumulative_var), group="") +
  ggtitle("percent variance explained, named, normalized root metabolites")

raw

leaf

met_amt.leaf_PCA <- met_amt %>%
  select(matches("_leaf_")) %>%
  prcomp(center = FALSE, scale. = FALSE) #already centered and scaled
names(met_per_mg.leaf_PCA)
[1] "sdev"     "rotation" "center"   "scale"    "x"       
tibble(variance=met_amt.leaf_PCA$sdev^2, PC=str_c("PC", 
                                                   str_pad(1:length(met_amt.leaf_PCA$sdev), width = 2, pad="0"))) %>%
  mutate(percent_var=100*variance/sum(variance),  
         cumulative_var=cumsum(percent_var)) %>%
  magrittr::extract(1:15,) %>%
  ggplot(aes(x=PC, y=percent_var)) +
  geom_col(fill="skyblue") + 
  geom_line(aes(y=cumulative_var), group="") +
  ggtitle("percent variance explained, named, raw leaf metabolites")

root

met_amt.root_PCA <- met_amt %>%
  select(matches("_root_")) %>%
  prcomp(center = FALSE, scale. = FALSE) #already centered and scaled
names(met_per_mg.root_PCA)
[1] "sdev"     "rotation" "center"   "scale"    "x"       
tibble(variance=met_amt.root_PCA$sdev^2, PC=str_c("PC", 
                                                   str_pad(1:length(met_amt.root_PCA$sdev), width = 2, pad="0"))) %>%
  mutate(percent_var=100*variance/sum(variance),  
         cumulative_var=cumsum(percent_var)) %>%
  magrittr::extract(1:15,) %>%
  ggplot(aes(x=PC, y=percent_var)) +
  geom_col(fill="skyblue") + 
  geom_line(aes(y=cumulative_var), group="") +
  ggtitle("percent variance explained, named, raw root metabolites")

now try these in a penalized regression

normalized

are the PCs normalized?

colMeans(met_amt.leaf_PCA$x) %>% round(3) #yes centered
 PC1  PC2  PC3  PC4  PC5  PC6  PC7  PC8  PC9 PC10 PC11 PC12 PC13 PC14 
   0    0    0    0    0    0    0    0    0    0    0    0    0    0 
PC15 PC16 PC17 PC18 PC19 PC20 PC21 PC22 PC23 PC24 
   0    0    0    0    0    0    0    0    0    0 
apply(met_amt.leaf_PCA$x, 2, sd) %>% round(2) #not scaled
  PC1   PC2   PC3   PC4   PC5   PC6   PC7   PC8   PC9  PC10  PC11  PC12 
11.89  9.03  6.73  6.47  5.80  5.78  5.37  5.08  4.58  4.40  4.36  4.35 
 PC13  PC14  PC15  PC16  PC17  PC18  PC19  PC20  PC21  PC22  PC23  PC24 
 4.03  3.83  3.76  3.64  3.55  3.46  3.30  3.19  3.14  3.06  0.00  0.00 

combine the leaf and root, and then scale them:

met_per_mg.leaf_PCs <- met_per_mg.leaf_PCA$x
colnames(met_per_mg.leaf_PCs) <- str_c("leaf_", colnames(met_per_mg.leaf_PCs))

met_per_mg.root_PCs <- met_per_mg.root_PCA$x
colnames(met_per_mg.root_PCs) <- str_c("root_", colnames(met_per_mg.root_PCs))

met_per_mg.PCs <- cbind(met_per_mg.leaf_PCs, met_per_mg.root_PCs) %>%
  scale()

met_amt.leaf_PCs <- met_amt.leaf_PCA$x
colnames(met_amt.leaf_PCs) <- str_c("leaf_", colnames(met_amt.leaf_PCs))

met_amt.root_PCs <- met_amt.root_PCA$x
colnames(met_amt.root_PCs) <- str_c("root_", colnames(met_amt.root_PCs))

met_amt.PCs <- cbind(met_amt.leaf_PCs, met_amt.root_PCs) %>%
  scale()

also combine the rotations

met_per_mg.leaf_rotation <- met_per_mg.leaf_PCA$rotation %>%
  as.data.frame() %>% 
  rename_with(~ str_c("leaf_", .x)) %>%
  rownames_to_column("metabolite")

met_per_mg.root_rotation <- met_per_mg.root_PCA$rotation %>%
  as.data.frame() %>% 
  rename_with(~ str_c("root_", .x)) %>%
  rownames_to_column("metabolite")

met_per_mg.PC_rotation <- full_join(met_per_mg.leaf_rotation, met_per_mg.root_rotation, by="metabolite")

met_amt.leaf_rotation <- met_amt.leaf_PCA$rotation %>% 
  as.data.frame() %>% 
  rename_with(~ str_c("leaf_", .x)) %>%
  rownames_to_column("metabolite")

met_amt.root_rotation <- met_amt.root_PCA$rotation %>%
  as.data.frame() %>% 
  rename_with(~ str_c("root_", .x)) %>%
  rownames_to_column("metabolite")

met_amt.PC_rotation <- full_join(met_amt.leaf_rotation, met_amt.root_rotation, by="metabolite")
met_per_mg_fit1LOO <- cv.glmnet(x=met_per_mg.PCs, y=leaflength$leaf_avg_std, nfolds = nrow(met_per_mg.PCs), alpha=1 )
Option grouped=FALSE enforced in cv.glmnet, since < 3 observations per fold
plot(met_per_mg_fit1LOO)

bestlam=met_per_mg_fit1LOO$lambda.1se

NEXT STEP: Do a K-fold CV, repeat many times and average. Might as well do alpha while we are at it. If we are doing alpha, then we need to manually create our own folds list for each run

normalized

multi CV

Fit 101 CVs for each of 11 alphas

set.seed(1245)

folds <- tibble(run=1:101) %>% 
  mutate(folds=map(run, ~ sample(rep(1:4,6))))

system.time (met_per_mg_multiCV <- expand_grid(run=1:100, alpha=round(seq(0,1,.1),1)) %>%
               left_join(folds, by="run") %>%
               mutate(fit=map2(folds, alpha, ~ cv.glmnet(x=met_per_mg.PCs, y=leaflength$leaf_avg_std, foldid = .x, alpha=.y
                                                         )))
             #, lambda=exp(seq(-5,0,length.out = 50)) )))
) #100 seconds
   user  system elapsed 
 42.060   1.285  44.588 
head(met_per_mg_multiCV)

for each fit, pull out the mean cv error, lambda, min lambda, and 1se lambda

met_per_mg_multiCV <- met_per_mg_multiCV %>%
  mutate(cvm=map(fit, magrittr::extract("cvm")),
         lambda=map(fit, magrittr::extract("lambda")),
         lambda.min=map_dbl(fit, magrittr::extract("lambda.min" )),
         lambda.1se=map_dbl(fit, magrittr::extract("lambda.1se")),
         nzero=map(fit, magrittr::extract("nzero"))
  )

head(met_per_mg_multiCV)

now calculate the mean and sem of cvm and min,1se labmdas. These need to be done separately because of the way the grouping works

met_per_mg_summary_cvm <- met_per_mg_multiCV %>% dplyr::select(-fit, -folds) %>% 
  unnest(c(cvm, lambda)) %>%
  group_by(alpha, lambda) %>%
  summarize(meancvm=mean(cvm), sem=sd(cvm)/sqrt(n()), high=meancvm+sem, low=meancvm-sem)
`summarise()` regrouping output by 'alpha' (override with `.groups` argument)
met_per_mg_summary_cvm
met_per_mg_summary_lambda <- met_per_mg_multiCV %>% dplyr::select(-fit, -folds, -cvm) %>% 
  group_by(alpha) %>%
  summarize(
    lambda.min.sd=sd(lambda.min), 
    lambda.min.mean=mean(lambda.min),
    #lambda.min.med=median(lambda.min), 
    lambda.min.high=lambda.min.mean+lambda.min.sd,
    #lambda.min.low=lambda.min.mean-lambda.min.sem,
    #lambda.1se.sem=sd(lambda.1se)/sqrt(n()), 
    lambda.1se.mean=mean(lambda.1se),
    #lambda.1se.med=median(lambda.1se), 
    #lambda.1se.high=lambda.1se+lambda.1se.sem,
    #lambda.1se.low=lambda.1se-lambda.1se.sem,
    nzero=nzero[1],
    lambda=lambda[1]
  )
`summarise()` ungrouping output (override with `.groups` argument)
met_per_mg_summary_lambda

plot it

met_per_mg_summary_cvm %>%
  #filter(alpha!=0) %>% # worse than everything else and throwing the plots off
  ggplot(aes(x=log(lambda), y= meancvm,  ymin=low, ymax=high)) +
  geom_ribbon(alpha=.25) +
  geom_line(aes(color=as.character(alpha))) +
  facet_wrap(~ as.character(alpha)) +
   coord_cartesian(xlim=(c(-5,0))) +
  geom_vline(aes(xintercept=log(lambda.min.mean)), alpha=.5, data=met_per_mg_summary_lambda) +
  geom_vline(aes(xintercept=log(lambda.min.high)), alpha=.5, data=met_per_mg_summary_lambda, color="blue") 

So overall these look more reasonable than the LOO plot.

Make a plot of MSE at minimum lambda for each alpha

met_per_mg_summary_cvm %>% 
  group_by(alpha) %>%
  filter(rank(meancvm, ties.method = "first")==1) %>%
  ggplot(aes(x=alpha,y=meancvm,ymin=low,ymax=high)) +
  geom_ribbon(color=NA, fill="gray80") +
  geom_line() +
  geom_point()

not a particular large difference there, aside from 0.1 and even then, not too much better.

Plot the number of nzero coefficients

met_per_mg_summary_lambda %>%
  unnest(c(lambda, nzero)) %>%
  group_by(alpha) %>%
  filter(abs(lambda.min.mean-lambda)==min(abs(lambda.min.mean-lambda))  ) %>%
  ungroup() %>%

ggplot(aes(x=as.character(alpha), y=nzero)) +
  geom_point() +
  ggtitle("Number of non-zero coefficents at minimum lambda") +
  ylim(0,24)

OK let’s do repeated test train starting from these CV lambdas

multi_tt <- function(lambda, alpha, n=10000, sample_size=24, train_size=20, x, y=leaflength$leaf_avg_std) {
  print(lambda)
  print(alpha)
tt <-
  tibble(run=1:n) %>%
  mutate(train=map(run, ~ sample(1:sample_size, train_size))) %>%
  mutate(fit=map(train, ~ glmnet(x=x[.,], y=y[.], lambda = lambda, alpha = alpha ))) %>%
  
  mutate(pred=map2(fit, train, ~ predict(.x, newx = x[-.y,]))) %>%
  mutate(cor=map2_dbl(pred, train, ~ cor(.x, y[-.y])  )) %>%
  mutate(MSE=map2_dbl(pred, train, ~ mean((y[-.y] - .x)^2))) %>%
  summarize(
    num_na=sum(is.na(cor)), 
    num_lt_0=sum(cor<=0, na.rm=TRUE),
    avg_cor=mean(cor, na.rm=TRUE),
    avg_MSE=mean(MSE))
tt
}

per_mg_fit_test_train <- met_per_mg_summary_lambda %>% 
  select(alpha, lambda.min.mean)

per_mg_fit_test_train <- met_per_mg_multiCV %>%
  filter(run==1) %>%
  select(alpha, fit) %>%
  right_join(per_mg_fit_test_train)
Joining, by = "alpha"
per_mg_fit_test_train <- per_mg_fit_test_train %>%
  mutate(pred_full=map2(fit, lambda.min.mean, ~ predict(.x, s=.y, newx=met_per_mg.PCs)),
         full_R=map_dbl(pred_full, ~ cor(.x, leaflength$leaf_avg_std)),
         full_MSE=map_dbl(pred_full, ~ mean((leaflength$leaf_avg_std-.x)^2))) %>%
  
  mutate(tt=map2(lambda.min.mean, alpha, ~ multi_tt(lambda=.x, alpha=.y, x=met_per_mg.PCs)))
[1] 63.76959
[1] 0
[1] 3.23426
[1] 0.1
[1] 1.96389
[1] 0.2
[1] 1.390028
[1] 0.3
[1] 1.057486
[1] 0.4
[1] 0.8344099
[1] 0.5
[1] 0.7009496
[1] 0.6
[1] 0.5978767
[1] 0.7
[1] 0.527807
[1] 0.8
[1] 0.4681233
[1] 0.9
[1] 0.4171747
[1] 1
(per_mg_fit_test_train <- per_mg_fit_test_train %>% unnest(tt))
per_mg_fit_test_train %>%
  ggplot(aes(x=alpha)) +
  geom_line(aes(y=avg_cor), color="red") +
  geom_point(aes(y=avg_cor), color="red") +
  geom_line(aes(y=avg_MSE), color="blue") +
  geom_point(aes(y=avg_MSE), color="blue")

all are similar from 0.2 onwards

look at fit:

alpha_per_mg <- .2

best_per_mg <- per_mg_fit_test_train %>% filter(alpha == alpha_per_mg) 
best_per_mg_fit <- best_per_mg$fit[[1]]
best_per_mg_lambda <- best_per_mg$lambda.min.mean

per_mg_coef.tb <- coef(best_per_mg_fit, s=best_per_mg_lambda) %>% 
  as.matrix() %>% as.data.frame() %>% 
  rownames_to_column(var="PC") %>%
  rename(beta=`1`)
  
per_mg_coef.tb %>% filter(beta!=0) %>% arrange(beta)
NA

pred and obs

plot(leaflength$leaf_avg_std, best_per_mg$pred_full[[1]])

cor.test(leaflength$leaf_avg_std, best_per_mg$pred_full[[1]]) #.57

    Pearson's product-moment correlation

data:  leaflength$leaf_avg_std and best_per_mg$pred_full[[1]]
t = 4.0007, df = 22, p-value = 0.0006021
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.3326307 0.8340147
sample estimates:
    cor 
0.64895 
best_per_mg$full_MSE
[1] 0.8477346

Percent variance explained

per_mg_vars <- per_mg_coef.tb %>% 
  filter(beta !=0, PC!="(Intercept)") %>%
  pull(PC) %>% c("leaf_avg_std", .)

per_mg_relimp <- leaflength %>% select(leaf_avg_std) %>% cbind(met_per_mg.PCs) %>% as.data.frame() %>% dplyr::select(all_of(per_mg_vars)) %>%
  calc.relimp() 

per_mg_coef.tb <- per_mg_relimp@lmg %>% as.matrix() %>% as.data.frame() %>%
  rownames_to_column("PC") %>%
  rename(PropVar_met_per_mg=V1) %>%
  full_join(per_mg_coef.tb) %>%
  arrange(desc(PropVar_met_per_mg))
Joining, by = "PC"
per_mg_coef.tb
NA

Checkout the rotations.

met_per_mg_rotation_out <- met_per_mg.PC_rotation %>% 
  pivot_longer(-metabolite, names_to="PC", values_to="loading") %>%
  filter(PC %in% filter(per_mg_coef.tb, beta!=0)$PC ) %>%
  group_by(PC) %>%
    filter(!str_detect(metabolite,".*(leaf|root)_[0-9]*$")) %>%
  filter(abs(loading) >= 0.05) %>%
  left_join(per_mg_coef.tb, by="PC") %>%
  arrange(desc(abs(beta)), desc(abs(loading))) %>%
  mutate(organ=ifelse(str_detect(metabolite, "_leaf_"), "leaf", "root"),
         transformation="normalized",
         metabolite=str_remove(metabolite, "met_per_mg_(root|leaf)_"),
         metabolite_effect_on_leaf=ifelse(beta*loading>0, "increase", "decrease"))
met_per_mg_rotation_out %>%  write_csv("../output/Leaf_associated_metabolites_normalized_NOBLANK.csv")
met_per_mg_rotation_out

non-normazlized

multi CV

Fit 101 CVs for each of 11 alphas

set.seed(1245)

folds <- tibble(run=1:101) %>% 
  mutate(folds=map(run, ~ sample(rep(1:6,4))))

system.time (met_amt_multiCV <- expand_grid(run=1:100, alpha=round(seq(0,1,.1),1)) %>%
               left_join(folds, by="run") %>%
               mutate(fit=map2(folds, alpha, ~ cv.glmnet(x=met_amt.PCs, y=leaflength$leaf_avg_std, foldid = .x, alpha=.y
                                                         )))
             #, lambda=exp(seq(-5,0,length.out = 50)) )))
) #100 seconds
   user  system elapsed 
 64.536   2.351  72.001 
head(met_amt_multiCV)

for each fit, pull out the mean cv error, lambda, min lambda, and 1se lambda

met_amt_multiCV <- met_amt_multiCV %>%
  mutate(cvm=map(fit, magrittr::extract("cvm")),
         lambda=map(fit, magrittr::extract("lambda")),
         lambda.min=map_dbl(fit, magrittr::extract("lambda.min" )),
         lambda.1se=map_dbl(fit, magrittr::extract("lambda.1se")),
         nzero=map(fit, magrittr::extract("nzero"))
  )

head(met_amt_multiCV)

now calculate the mean and sem of cvm and min,1se labmdas. These need to be done separately because of the way the grouping works

met_amt_summary_cvm <- met_amt_multiCV %>% dplyr::select(-fit, -folds) %>% 
  unnest(c(cvm, lambda)) %>%
  group_by(alpha, lambda) %>%
  summarize(meancvm=mean(cvm), sem=sd(cvm)/sqrt(n()), high=meancvm+sem, low=meancvm-sem)
`summarise()` regrouping output by 'alpha' (override with `.groups` argument)
met_amt_summary_cvm
met_amt_summary_lambda <- met_amt_multiCV %>% dplyr::select(-fit, -folds, -cvm) %>% 
  group_by(alpha) %>%
  summarize(
    lambda.min.sd=sd(lambda.min), 
    lambda.min.mean=mean(lambda.min),
    #lambda.min.med=median(lambda.min), 
    lambda.min.high=lambda.min.mean+lambda.min.sd,
    #lambda.min.low=lambda.min.mean-lambda.min.sem,
    #lambda.1se.sem=sd(lambda.1se)/sqrt(n()), 
    lambda.1se.mean=mean(lambda.1se),
    #lambda.1se.med=median(lambda.1se), 
    #lambda.1se.high=lambda.1se+lambda.1se.sem,
    #lambda.1se.low=lambda.1se-lambda.1se.sem,
    nzero=nzero[1],
    lambda=lambda[1]
  )
`summarise()` ungrouping output (override with `.groups` argument)
met_amt_summary_lambda

plot it

met_amt_summary_cvm %>%
  #filter(alpha!=0) %>% # worse than everything else and throwing the plots off
  ggplot(aes(x=log(lambda), y= meancvm,  ymin=low, ymax=high)) +
  geom_ribbon(alpha=.25) +
  geom_line(aes(color=as.character(alpha))) +
  facet_wrap(~ as.character(alpha)) +
   coord_cartesian(xlim=(c(-5,0))) +
  geom_vline(aes(xintercept=log(lambda.min.mean)), alpha=.5, data=met_amt_summary_lambda) +
  geom_vline(aes(xintercept=log(lambda.min.high)), alpha=.5, data=met_amt_summary_lambda, color="blue") 

Make a plot of MSE at minimum lambda for each alpha

met_amt_summary_cvm %>% 
  group_by(alpha) %>%
  filter(rank(meancvm, ties.method = "first")==1) %>%
  ggplot(aes(x=alpha,y=meancvm,ymin=low,ymax=high)) +
  geom_ribbon(color=NA, fill="gray80") +
  geom_line() +
  geom_point()

decreasing throughout, although differences are pretty small

Plot the number of nzero coefficients

met_amt_summary_lambda %>%
  unnest(c(lambda, nzero)) %>%
  group_by(alpha) %>%
  filter(abs(lambda.min.mean-lambda)==min(abs(lambda.min.mean-lambda))  ) %>%
  ungroup() %>%

ggplot(aes(x=as.character(alpha), y=nzero)) +
  geom_point() +
  ggtitle("Number of non-zero coefficents at minimum lambda") +
  ylim(0,36)

OK let’s do repeated test train starting from these CV lambdas

multi_tt <- function(lambda, alpha, n=10000, sample_size=24, train_size=20, x, y=leaflength$leaf_avg_std) {
  print(lambda)
  print(alpha)
tt <-
  tibble(run=1:n) %>%
  mutate(train=map(run, ~ sample(1:sample_size, train_size))) %>%
  mutate(fit=map(train, ~ glmnet(x=x[.,], y=y[.], lambda = lambda, alpha = alpha ))) %>%
  
  mutate(pred=map2(fit, train, ~ predict(.x, newx = x[-.y,]))) %>%
  mutate(cor=map2_dbl(pred, train, ~ cor(.x, y[-.y])  )) %>%
  mutate(MSE=map2_dbl(pred, train, ~ mean((y[-.y] - .x)^2))) %>%
  summarize(
    num_na=sum(is.na(cor)), 
    num_lt_0=sum(cor<=0, na.rm=TRUE),
    avg_cor=mean(cor, na.rm=TRUE),
    avg_MSE=mean(MSE))
tt
}

amt_fit_test_train <- met_amt_summary_lambda %>% 
  select(alpha, lambda.min.mean)

amt_fit_test_train <- met_amt_multiCV %>%
  filter(run==1) %>%
  select(alpha, fit) %>%
  right_join(amt_fit_test_train)
Joining, by = "alpha"
amt_fit_test_train <- amt_fit_test_train %>%
  mutate(pred_full=map2(fit, lambda.min.mean, ~ predict(.x, s=.y, newx=met_amt.PCs)),
         full_R=map_dbl(pred_full, ~ cor(.x, leaflength$leaf_avg_std)),
         full_MSE=map_dbl(pred_full, ~ mean((leaflength$leaf_avg_std-.x)^2))) %>%
  
  mutate(tt=map2(lambda.min.mean, alpha, ~ multi_tt(lambda=.x, alpha=.y, x=met_amt.PCs)))
[1] 129.3889
[1] 0
[1] 2.394178
[1] 0.1
[1] 1.622926
[1] 0.2
[1] 1.129161
[1] 0.3
[1] 0.8782994
[1] 0.4
[1] 0.709034
[1] 0.5
[1] 0.5979773
[1] 0.6
[1] 0.5183774
[1] 0.7
[1] 0.4567253
[1] 0.8
[1] 0.4084716
[1] 0.9
[1] 0.3686524
[1] 1
(amt_fit_test_train <- amt_fit_test_train %>% unnest(tt))
amt_fit_test_train %>%
  ggplot(aes(x=alpha)) +
  geom_line(aes(y=avg_cor), color="red") +
  geom_point(aes(y=avg_cor), color="red") +
  geom_line(aes(y=avg_MSE), color="blue") +
  geom_point(aes(y=avg_MSE), color="blue")

not a big difference after 0.2

look at fit:

alpha_amt <- .7

best_amt <- amt_fit_test_train %>% filter(alpha == alpha_amt) 
best_amt_fit <- best_amt$fit[[1]]
best_amt_lambda <- best_amt$lambda.min.mean

amt_coef.tb <- coef(best_amt_fit, s=best_amt_lambda) %>% 
  as.matrix() %>% as.data.frame() %>% 
  rownames_to_column(var="PC") %>%
  rename(beta=`1`)
  
amt_coef.tb %>% filter(beta!=0) %>% arrange(beta)
NA

pred and obs

plot(leaflength$leaf_avg_std, best_amt$pred_full[[1]])

cor.test(leaflength$leaf_avg_std, best_amt$pred_full[[1]]) #.64

    Pearson's product-moment correlation

data:  leaflength$leaf_avg_std and best_amt$pred_full[[1]]
t = 3.9359, df = 22, p-value = 0.000705
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.3232239 0.8307777
sample estimates:
      cor 
0.6428067 
best_amt$full_MSE
[1] 0.6804825

Percent variance explained

amt_vars <- amt_coef.tb %>% 
  filter(beta !=0, PC!="(Intercept)") %>%
  pull(PC) %>% c("leaf_avg_std", .)

# because there is only one explanatory variable, we don't need to (and can't) use relimp
amt_relimp <- leaflength %>% select(leaf_avg_std) %>% cbind(met_amt.PCs) %>% as.data.frame() %>% dplyr::select(all_of(amt_vars)) %>% cor() %>% magrittr::extract(1,2) %>% magrittr::multiply_by(.,.)
amt_relimp <- tibble(PC=str_subset(amt_vars, "PC"),
                     PropVar_met_amt=amt_relimp)


amt_coef.tb <- amt_relimp %>% 
  full_join(amt_coef.tb) %>%
  arrange(desc(PropVar_met_amt))
Joining, by = "PC"
amt_coef.tb
NA

Checkout the rotations.

met_amt_rotation_out <- met_amt.PC_rotation %>% 
  pivot_longer(-metabolite, names_to="PC", values_to="loading") %>%
  filter(PC %in% filter(amt_coef.tb, beta!=0)$PC ) %>%
  group_by(PC) %>%
    filter(!str_detect(metabolite,".*(leaf|root)_[0-9]*$")) %>%
  filter(abs(loading) >= 0.05) %>%
  left_join(amt_coef.tb, by="PC") %>%
  arrange(desc(abs(beta)), desc(abs(loading))) %>%
  mutate(organ=ifelse(str_detect(metabolite, "_leaf_"), "leaf", "root"),
         transformation="raw",
         metabolite=str_remove(metabolite, "met_amt_(root|leaf)_"),
         metabolite_effect_on_leaf=ifelse(beta*loading>0, "increase", "decrease"))
met_amt_rotation_out %>%  write_csv("../output/Leaf_associated_metabolites_raw_NOBLANK.csv")

met_amt_rotation_out
LS0tCnRpdGxlOiAiTWFjaGluZSBsZWFybmluZyBmb3IgbWV0YWJvbGl0ZXMsIG5vIGJsYW5rLCBhbGwgbWV0YWJvbGl0ZXMiCmF1dGhvcjogIkp1bGluIE1hbG9vZiIKZGF0ZTogIjEvMjgvMjAyMSIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKYGBge3Igc2V0dXAsIGluY2x1ZGU9RkFMU0V9CmtuaXRyOjpvcHRzX2NodW5rJHNldChlY2hvID0gVFJVRSkKYGBgCgpHb2FsOiBjYW4gbWFjaGluZSBsZWFybmluZyBtZXRob2RzIGhlbHAgdXMgdG8gYXNzb2NpYXRlIG1ldGFib2xpdGVzIHdpdGggbGVhZiBsZW5ndGg/CgpQcmV2aW91c2x5IChzY3JpcHQgMTFiMikgSSBmaWx0ZXJlZCBvdXQgdW5uYW1lZCBtZXRhYm9saXRlcy4gIEhlcmUgSSBrZWVwIHRoZW0gYWxsLiAgYnV0IHJlbW92ZSBibGFuayBzYW1wbGVzCgpgYGB7cn0KbGlicmFyeShnbG1uZXQpCmxpYnJhcnkocmVsYWltcG8pCmxpYnJhcnkodGlkeXZlcnNlKQpgYGAKCmdldCBsZWFmbGVuZ3RoIGRhdGEKYGBge3J9CmxlYWZsZW5ndGggPC0gcmVhZF9jc3YoIi4uLy4uL3BsYW50L291dHB1dC9sZWFmX2xlbmd0aHNfbWV0YWJvbGl0ZS5jc3YiKSAlPiUKICBtdXRhdGUocG90PXN0cl9wYWQocG90LCB3aWR0aD0zLCBwYWQ9IjAiKSwKICAgICAgICAgc2FtcGxlSUQ9c3RyX2MoInd5byIsIGdlbm90eXBlLCBwb3QsIHNlcD0iXyIpLCAKICAgICAgICAgZ3JvdXA9c3RyX2Moc29pbCxnZW5vdHlwZSx0cnQsc2VwPSJfIikpICU+JQogIHNlbGVjdChzYW1wbGVJRCwgZ3JvdXAsIGxlYWZfYXZnX3N0ZCkgICU+JQogIGZpbHRlcighc3RyX2RldGVjdChncm91cCwgIkJMQU5LIikpCmxlbmd0aCh1bmlxdWUobGVhZmxlbmd0aCRncm91cCkpCmxlYWZsZW5ndGggJT4lIGFycmFuZ2Uoc2FtcGxlSUQpCmBgYAoKZ2V0IGFuZCB3cmFuZ2xlIG1ldGFib2xpdGUgZGF0YQpgYGB7cn0KbWV0X3JhdyA8LXJlYWRfY3N2KCIuLi9pbnB1dC9tZXRhYm9saXRlc19zZXQxLmNzdiIpCm1ldCA8LSBtZXRfcmF3ICU+JSAKICBtdXRhdGUocG90PXN0cl9wYWQocG90LCB3aWR0aCA9IDMsIHBhZCA9ICIwIikpICU+JQogIG11dGF0ZShzYW1wbGVJRD1zdHJfYygid3lvIiwgZ2Vub3R5cGUsIHBvdCwgc2VwPSJfIikpICU+JQogIGZpbHRlcihzb2lsIT0iQkxBTksiKSAlPiUKICBzZWxlY3Qoc2FtcGxlSUQsIGdlbm90eXBlLCB0aXNzdWUsIHNhbXBsZV9tYXNzID0gYHNhbXBsZV9tYXNzIG1nYCwgIXN1Ym1pc3Npb25fbnVtYmVyOmNvbmNhdGVuYXRlKSAlPiUKICBwaXZvdF9sb25nZXIoIXNhbXBsZUlEOnNhbXBsZV9tYXNzLCBuYW1lc190byA9ICJtZXRhYm9saXRlIiwgdmFsdWVzX3RvID0gIm1ldF9hbW91bnQiKSAlPiUKICAKICAjYWRqdXN0IGJ5IHNhbXBsZSBtYXNzCiAgbXV0YXRlKG1ldF9wZXJfbWc9bWV0X2Ftb3VudC9zYW1wbGVfbWFzcykgJT4lCiAgCiAgI3NjYWxlIGFuZCBjZW50ZXIKICBncm91cF9ieShtZXRhYm9saXRlLCBnZW5vdHlwZSwgdGlzc3VlKSAlPiUKICBtdXRhdGUobWV0X3Blcl9tZz1zY2FsZShtZXRfcGVyX21nKSwKICAgICAgICAgbWV0X2FtdD1zY2FsZShtZXRfYW1vdW50KQogICkgJT4lIAogIHBpdm90X3dpZGVyKGlkX2NvbHMgPSBzYW1wbGVJRCwgCiAgICAgICAgICAgICAgbmFtZXNfZnJvbSA9IGModGlzc3VlLCBtZXRhYm9saXRlKSwgCiAgICAgICAgICAgICAgdmFsdWVzX2Zyb20gPSBzdGFydHNfd2l0aCgibWV0XyIpLAogICAgICAgICAgICAgIG5hbWVzX3NlcCA9ICJfIikKCm1ldApgYGAKCnNwbGl0IHRoaXMgaW50byB0d28gZGF0YSBmcmFtZXMsIG9uZSBub3JtYWxpemVkIGJ5IHRpc3N1ZSBhbW91bnQgYW5kIG9uZSBub3QuCmBgYHtyfQptZXRfcGVyX21nIDwtIG1ldCAlPiUgc2VsZWN0KHNhbXBsZUlELCAgc3RhcnRzX3dpdGgoIm1ldF9wZXJfbWciKSkgJT4lCiAgYXMuZGF0YS5mcmFtZSgpICU+JSBjb2x1bW5fdG9fcm93bmFtZXMoInNhbXBsZUlEIikKbWV0X2FtdCA8LSBtZXQgJT4lIHNlbGVjdChzYW1wbGVJRCwgIHN0YXJ0c193aXRoKCJtZXRfYW10IikpICU+JQogIGFzLmRhdGEuZnJhbWUoKSAlPiUgY29sdW1uX3RvX3Jvd25hbWVzKCJzYW1wbGVJRCIpCmBgYAoKZ2V0IGxlYWYgZGF0YSBvcmRlciB0byBtYXRjaAoKYGBge3J9CmxlYWZsZW5ndGggPC0gbGVhZmxlbmd0aFttYXRjaChtZXQkc2FtcGxlSUQsIGxlYWZsZW5ndGgkc2FtcGxlSUQpLF0KbGVhZmxlbmd0aApgYGAKCgojIyBDYWxjIFBDQXM6CgojIyMgbm9ybWFsaXplZAoKIyMjIyBsZWFmCgpgYGB7ciwgZmlnLmFzcD0xfQptZXRfcGVyX21nLmxlYWZfUENBIDwtIG1ldF9wZXJfbWcgJT4lIAogIHNlbGVjdChtYXRjaGVzKCJfbGVhZl8iKSkgJT4lCiAgcHJjb21wKGNlbnRlciA9IEZBTFNFLCBzY2FsZS4gPSBGQUxTRSkgI2FscmVhZHkgY2VudGVyZWQgYW5kIHNjYWxlZApuYW1lcyhtZXRfcGVyX21nLmxlYWZfUENBKQp0aWJibGUodmFyaWFuY2U9bWV0X3Blcl9tZy5sZWFmX1BDQSRzZGV2XjIsIFBDPXN0cl9jKCJQQyIsIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBzdHJfcGFkKDE6bGVuZ3RoKG1ldF9wZXJfbWcubGVhZl9QQ0Ekc2RldiksIHdpZHRoID0gMiwgcGFkPSIwIikpKSAlPiUKICBtdXRhdGUocGVyY2VudF92YXI9MTAwKnZhcmlhbmNlL3N1bSh2YXJpYW5jZSksICAKICAgICAgICAgY3VtdWxhdGl2ZV92YXI9Y3Vtc3VtKHBlcmNlbnRfdmFyKSkgJT4lCiAgbWFncml0dHI6OmV4dHJhY3QoMToxNSwpICU+JQogIGdncGxvdChhZXMoeD1QQywgeT1wZXJjZW50X3ZhcikpICsKICBnZW9tX2NvbChmaWxsPSJza3libHVlIikgKyAKICBnZW9tX2xpbmUoYWVzKHk9Y3VtdWxhdGl2ZV92YXIpLCBncm91cD0iIikgKwogIGdndGl0bGUoInBlcmNlbnQgdmFyaWFuY2UgZXhwbGFpbmVkLCBuYW1lZCwgbm9ybWFsaXplZCBsZWFmIG1ldGFib2xpdGVzIikKYGBgCgojIyMjIHJvb3QKCmBgYHtyLCBmaWcuYXNwPTF9Cm1ldF9wZXJfbWcucm9vdF9QQ0EgPC0gbWV0X3Blcl9tZyAlPiUgCiAgc2VsZWN0KG1hdGNoZXMoIl9yb290XyIpKSAlPiUKICBwcmNvbXAoY2VudGVyID0gRkFMU0UsIHNjYWxlLiA9IEZBTFNFKSAjYWxyZWFkeSBjZW50ZXJlZCBhbmQgc2NhbGVkCm5hbWVzKG1ldF9wZXJfbWcucm9vdF9QQ0EpCnRpYmJsZSh2YXJpYW5jZT1tZXRfcGVyX21nLnJvb3RfUENBJHNkZXZeMiwgUEM9c3RyX2MoIlBDIiwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHN0cl9wYWQoMTpsZW5ndGgobWV0X3Blcl9tZy5yb290X1BDQSRzZGV2KSwgd2lkdGggPSAyLCBwYWQ9IjAiKSkpICU+JQogIG11dGF0ZShwZXJjZW50X3Zhcj0xMDAqdmFyaWFuY2Uvc3VtKHZhcmlhbmNlKSwgIAogICAgICAgICBjdW11bGF0aXZlX3Zhcj1jdW1zdW0ocGVyY2VudF92YXIpKSAlPiUKICBtYWdyaXR0cjo6ZXh0cmFjdCgxOjE1LCkgJT4lCiAgZ2dwbG90KGFlcyh4PVBDLCB5PXBlcmNlbnRfdmFyKSkgKwogIGdlb21fY29sKGZpbGw9InNreWJsdWUiKSArIAogIGdlb21fbGluZShhZXMoeT1jdW11bGF0aXZlX3ZhciksIGdyb3VwPSIiKSArCiAgZ2d0aXRsZSgicGVyY2VudCB2YXJpYW5jZSBleHBsYWluZWQsIG5hbWVkLCBub3JtYWxpemVkIHJvb3QgbWV0YWJvbGl0ZXMiKQpgYGAKCiMjIyByYXcKCiMjIyMgbGVhZgpgYGB7ciwgZmlnLmFzcD0xfQptZXRfYW10LmxlYWZfUENBIDwtIG1ldF9hbXQgJT4lCiAgc2VsZWN0KG1hdGNoZXMoIl9sZWFmXyIpKSAlPiUKICBwcmNvbXAoY2VudGVyID0gRkFMU0UsIHNjYWxlLiA9IEZBTFNFKSAjYWxyZWFkeSBjZW50ZXJlZCBhbmQgc2NhbGVkCm5hbWVzKG1ldF9wZXJfbWcubGVhZl9QQ0EpCnRpYmJsZSh2YXJpYW5jZT1tZXRfYW10LmxlYWZfUENBJHNkZXZeMiwgUEM9c3RyX2MoIlBDIiwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHN0cl9wYWQoMTpsZW5ndGgobWV0X2FtdC5sZWFmX1BDQSRzZGV2KSwgd2lkdGggPSAyLCBwYWQ9IjAiKSkpICU+JQogIG11dGF0ZShwZXJjZW50X3Zhcj0xMDAqdmFyaWFuY2Uvc3VtKHZhcmlhbmNlKSwgIAogICAgICAgICBjdW11bGF0aXZlX3Zhcj1jdW1zdW0ocGVyY2VudF92YXIpKSAlPiUKICBtYWdyaXR0cjo6ZXh0cmFjdCgxOjE1LCkgJT4lCiAgZ2dwbG90KGFlcyh4PVBDLCB5PXBlcmNlbnRfdmFyKSkgKwogIGdlb21fY29sKGZpbGw9InNreWJsdWUiKSArIAogIGdlb21fbGluZShhZXMoeT1jdW11bGF0aXZlX3ZhciksIGdyb3VwPSIiKSArCiAgZ2d0aXRsZSgicGVyY2VudCB2YXJpYW5jZSBleHBsYWluZWQsIG5hbWVkLCByYXcgbGVhZiBtZXRhYm9saXRlcyIpCmBgYAoKIyMjIyByb290CmBgYHtyLCBmaWcuYXNwPTF9Cm1ldF9hbXQucm9vdF9QQ0EgPC0gbWV0X2FtdCAlPiUKICBzZWxlY3QobWF0Y2hlcygiX3Jvb3RfIikpICU+JQogIHByY29tcChjZW50ZXIgPSBGQUxTRSwgc2NhbGUuID0gRkFMU0UpICNhbHJlYWR5IGNlbnRlcmVkIGFuZCBzY2FsZWQKbmFtZXMobWV0X3Blcl9tZy5yb290X1BDQSkKdGliYmxlKHZhcmlhbmNlPW1ldF9hbXQucm9vdF9QQ0Ekc2Rldl4yLCBQQz1zdHJfYygiUEMiLCAKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgc3RyX3BhZCgxOmxlbmd0aChtZXRfYW10LnJvb3RfUENBJHNkZXYpLCB3aWR0aCA9IDIsIHBhZD0iMCIpKSkgJT4lCiAgbXV0YXRlKHBlcmNlbnRfdmFyPTEwMCp2YXJpYW5jZS9zdW0odmFyaWFuY2UpLCAgCiAgICAgICAgIGN1bXVsYXRpdmVfdmFyPWN1bXN1bShwZXJjZW50X3ZhcikpICU+JQogIG1hZ3JpdHRyOjpleHRyYWN0KDE6MTUsKSAlPiUKICBnZ3Bsb3QoYWVzKHg9UEMsIHk9cGVyY2VudF92YXIpKSArCiAgZ2VvbV9jb2woZmlsbD0ic2t5Ymx1ZSIpICsgCiAgZ2VvbV9saW5lKGFlcyh5PWN1bXVsYXRpdmVfdmFyKSwgZ3JvdXA9IiIpICsKICBnZ3RpdGxlKCJwZXJjZW50IHZhcmlhbmNlIGV4cGxhaW5lZCwgbmFtZWQsIHJhdyByb290IG1ldGFib2xpdGVzIikKYGBgCgojIyBub3cgdHJ5IHRoZXNlIGluIGEgcGVuYWxpemVkIHJlZ3Jlc3Npb24KCiMjIyBub3JtYWxpemVkCgphcmUgdGhlIFBDcyBub3JtYWxpemVkPwpgYGB7cn0KY29sTWVhbnMobWV0X2FtdC5sZWFmX1BDQSR4KSAlPiUgcm91bmQoMykgI3llcyBjZW50ZXJlZAphcHBseShtZXRfYW10LmxlYWZfUENBJHgsIDIsIHNkKSAlPiUgcm91bmQoMikgI25vdCBzY2FsZWQKYGBgCgpjb21iaW5lIHRoZSBsZWFmIGFuZCByb290LCBhbmQgdGhlbiBzY2FsZSB0aGVtOgpgYGB7cn0KbWV0X3Blcl9tZy5sZWFmX1BDcyA8LSBtZXRfcGVyX21nLmxlYWZfUENBJHgKY29sbmFtZXMobWV0X3Blcl9tZy5sZWFmX1BDcykgPC0gc3RyX2MoImxlYWZfIiwgY29sbmFtZXMobWV0X3Blcl9tZy5sZWFmX1BDcykpCgptZXRfcGVyX21nLnJvb3RfUENzIDwtIG1ldF9wZXJfbWcucm9vdF9QQ0EkeApjb2xuYW1lcyhtZXRfcGVyX21nLnJvb3RfUENzKSA8LSBzdHJfYygicm9vdF8iLCBjb2xuYW1lcyhtZXRfcGVyX21nLnJvb3RfUENzKSkKCm1ldF9wZXJfbWcuUENzIDwtIGNiaW5kKG1ldF9wZXJfbWcubGVhZl9QQ3MsIG1ldF9wZXJfbWcucm9vdF9QQ3MpICU+JQogIHNjYWxlKCkKCm1ldF9hbXQubGVhZl9QQ3MgPC0gbWV0X2FtdC5sZWFmX1BDQSR4CmNvbG5hbWVzKG1ldF9hbXQubGVhZl9QQ3MpIDwtIHN0cl9jKCJsZWFmXyIsIGNvbG5hbWVzKG1ldF9hbXQubGVhZl9QQ3MpKQoKbWV0X2FtdC5yb290X1BDcyA8LSBtZXRfYW10LnJvb3RfUENBJHgKY29sbmFtZXMobWV0X2FtdC5yb290X1BDcykgPC0gc3RyX2MoInJvb3RfIiwgY29sbmFtZXMobWV0X2FtdC5yb290X1BDcykpCgptZXRfYW10LlBDcyA8LSBjYmluZChtZXRfYW10LmxlYWZfUENzLCBtZXRfYW10LnJvb3RfUENzKSAlPiUKICBzY2FsZSgpCmBgYAoKYWxzbyBjb21iaW5lIHRoZSByb3RhdGlvbnMKYGBge3J9Cm1ldF9wZXJfbWcubGVhZl9yb3RhdGlvbiA8LSBtZXRfcGVyX21nLmxlYWZfUENBJHJvdGF0aW9uICU+JQogIGFzLmRhdGEuZnJhbWUoKSAlPiUgCiAgcmVuYW1lX3dpdGgofiBzdHJfYygibGVhZl8iLCAueCkpICU+JQogIHJvd25hbWVzX3RvX2NvbHVtbigibWV0YWJvbGl0ZSIpCgptZXRfcGVyX21nLnJvb3Rfcm90YXRpb24gPC0gbWV0X3Blcl9tZy5yb290X1BDQSRyb3RhdGlvbiAlPiUKICBhcy5kYXRhLmZyYW1lKCkgJT4lIAogIHJlbmFtZV93aXRoKH4gc3RyX2MoInJvb3RfIiwgLngpKSAlPiUKICByb3duYW1lc190b19jb2x1bW4oIm1ldGFib2xpdGUiKQoKbWV0X3Blcl9tZy5QQ19yb3RhdGlvbiA8LSBmdWxsX2pvaW4obWV0X3Blcl9tZy5sZWFmX3JvdGF0aW9uLCBtZXRfcGVyX21nLnJvb3Rfcm90YXRpb24sIGJ5PSJtZXRhYm9saXRlIikKCm1ldF9hbXQubGVhZl9yb3RhdGlvbiA8LSBtZXRfYW10LmxlYWZfUENBJHJvdGF0aW9uICU+JSAKICBhcy5kYXRhLmZyYW1lKCkgJT4lIAogIHJlbmFtZV93aXRoKH4gc3RyX2MoImxlYWZfIiwgLngpKSAlPiUKICByb3duYW1lc190b19jb2x1bW4oIm1ldGFib2xpdGUiKQoKbWV0X2FtdC5yb290X3JvdGF0aW9uIDwtIG1ldF9hbXQucm9vdF9QQ0Ekcm90YXRpb24gJT4lCiAgYXMuZGF0YS5mcmFtZSgpICU+JSAKICByZW5hbWVfd2l0aCh+IHN0cl9jKCJyb290XyIsIC54KSkgJT4lCiAgcm93bmFtZXNfdG9fY29sdW1uKCJtZXRhYm9saXRlIikKCm1ldF9hbXQuUENfcm90YXRpb24gPC0gZnVsbF9qb2luKG1ldF9hbXQubGVhZl9yb3RhdGlvbiwgbWV0X2FtdC5yb290X3JvdGF0aW9uLCBieT0ibWV0YWJvbGl0ZSIpCgpgYGAKCgoKYGBge3J9Cm1ldF9wZXJfbWdfZml0MUxPTyA8LSBjdi5nbG1uZXQoeD1tZXRfcGVyX21nLlBDcywgeT1sZWFmbGVuZ3RoJGxlYWZfYXZnX3N0ZCwgbmZvbGRzID0gbnJvdyhtZXRfcGVyX21nLlBDcyksIGFscGhhPTEgKQpwbG90KG1ldF9wZXJfbWdfZml0MUxPTykKYmVzdGxhbT1tZXRfcGVyX21nX2ZpdDFMT08kbGFtYmRhLjFzZQpgYGAKCiAgCgpORVhUIFNURVA6IERvIGEgSy1mb2xkIENWLCByZXBlYXQgbWFueSB0aW1lcyBhbmQgYXZlcmFnZS4gIE1pZ2h0IGFzIHdlbGwgZG8gYWxwaGEgd2hpbGUgd2UgYXJlIGF0IGl0LiAgSWYgd2UgYXJlIGRvaW5nIGFscGhhLCB0aGVuIHdlIG5lZWQgdG8gbWFudWFsbHkgY3JlYXRlIG91ciBvd24gZm9sZHMgbGlzdCBmb3IgZWFjaCBydW4KCiMgbm9ybWFsaXplZAoKIyMgbXVsdGkgQ1YKCkZpdCAxMDEgQ1ZzIGZvciBlYWNoIG9mIDExIGFscGhhcwpgYGB7cn0Kc2V0LnNlZWQoMTI0NSkKCmZvbGRzIDwtIHRpYmJsZShydW49MToxMDEpICU+JSAKICBtdXRhdGUoZm9sZHM9bWFwKHJ1biwgfiBzYW1wbGUocmVwKDE6NCw2KSkpKQoKc3lzdGVtLnRpbWUgKG1ldF9wZXJfbWdfbXVsdGlDViA8LSBleHBhbmRfZ3JpZChydW49MToxMDAsIGFscGhhPXJvdW5kKHNlcSgwLDEsLjEpLDEpKSAlPiUKICAgICAgICAgICAgICAgbGVmdF9qb2luKGZvbGRzLCBieT0icnVuIikgJT4lCiAgICAgICAgICAgICAgIG11dGF0ZShmaXQ9bWFwMihmb2xkcywgYWxwaGEsIH4gY3YuZ2xtbmV0KHg9bWV0X3Blcl9tZy5QQ3MsIHk9bGVhZmxlbmd0aCRsZWFmX2F2Z19zdGQsIGZvbGRpZCA9IC54LCBhbHBoYT0ueQogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICApKSkKICAgICAgICAgICAgICMsIGxhbWJkYT1leHAoc2VxKC01LDAsbGVuZ3RoLm91dCA9IDUwKSkgKSkpCikgIzEwMCBzZWNvbmRzCgpoZWFkKG1ldF9wZXJfbWdfbXVsdGlDVikKYGBgCgpmb3IgZWFjaCBmaXQsIHB1bGwgb3V0IHRoZSBtZWFuIGN2IGVycm9yLCBsYW1iZGEsIG1pbiBsYW1iZGEsIGFuZCAxc2UgbGFtYmRhIApgYGB7cn0KbWV0X3Blcl9tZ19tdWx0aUNWIDwtIG1ldF9wZXJfbWdfbXVsdGlDViAlPiUKICBtdXRhdGUoY3ZtPW1hcChmaXQsIG1hZ3JpdHRyOjpleHRyYWN0KCJjdm0iKSksCiAgICAgICAgIGxhbWJkYT1tYXAoZml0LCBtYWdyaXR0cjo6ZXh0cmFjdCgibGFtYmRhIikpLAogICAgICAgICBsYW1iZGEubWluPW1hcF9kYmwoZml0LCBtYWdyaXR0cjo6ZXh0cmFjdCgibGFtYmRhLm1pbiIgKSksCiAgICAgICAgIGxhbWJkYS4xc2U9bWFwX2RibChmaXQsIG1hZ3JpdHRyOjpleHRyYWN0KCJsYW1iZGEuMXNlIikpLAogICAgICAgICBuemVybz1tYXAoZml0LCBtYWdyaXR0cjo6ZXh0cmFjdCgibnplcm8iKSkKICApCgpoZWFkKG1ldF9wZXJfbWdfbXVsdGlDVikKYGBgCgoKbm93IGNhbGN1bGF0ZSB0aGUgbWVhbiBhbmQgc2VtIG9mIGN2bSBhbmQgbWluLDFzZSBsYWJtZGFzLiAgVGhlc2UgbmVlZCB0byBiZSBkb25lIHNlcGFyYXRlbHkgYmVjYXVzZSBvZiB0aGUgd2F5IHRoZSBncm91cGluZyB3b3JrcwpgYGB7cn0KbWV0X3Blcl9tZ19zdW1tYXJ5X2N2bSA8LSBtZXRfcGVyX21nX211bHRpQ1YgJT4lIGRwbHlyOjpzZWxlY3QoLWZpdCwgLWZvbGRzKSAlPiUgCiAgdW5uZXN0KGMoY3ZtLCBsYW1iZGEpKSAlPiUKICBncm91cF9ieShhbHBoYSwgbGFtYmRhKSAlPiUKICBzdW1tYXJpemUobWVhbmN2bT1tZWFuKGN2bSksIHNlbT1zZChjdm0pL3NxcnQobigpKSwgaGlnaD1tZWFuY3ZtK3NlbSwgbG93PW1lYW5jdm0tc2VtKQoKbWV0X3Blcl9tZ19zdW1tYXJ5X2N2bQpgYGAKCmBgYHtyfQptZXRfcGVyX21nX3N1bW1hcnlfbGFtYmRhIDwtIG1ldF9wZXJfbWdfbXVsdGlDViAlPiUgZHBseXI6OnNlbGVjdCgtZml0LCAtZm9sZHMsIC1jdm0pICU+JSAKICBncm91cF9ieShhbHBoYSkgJT4lCiAgc3VtbWFyaXplKAogICAgbGFtYmRhLm1pbi5zZD1zZChsYW1iZGEubWluKSwgCiAgICBsYW1iZGEubWluLm1lYW49bWVhbihsYW1iZGEubWluKSwKICAgICNsYW1iZGEubWluLm1lZD1tZWRpYW4obGFtYmRhLm1pbiksIAogICAgbGFtYmRhLm1pbi5oaWdoPWxhbWJkYS5taW4ubWVhbitsYW1iZGEubWluLnNkLAogICAgI2xhbWJkYS5taW4ubG93PWxhbWJkYS5taW4ubWVhbi1sYW1iZGEubWluLnNlbSwKICAgICNsYW1iZGEuMXNlLnNlbT1zZChsYW1iZGEuMXNlKS9zcXJ0KG4oKSksIAogICAgbGFtYmRhLjFzZS5tZWFuPW1lYW4obGFtYmRhLjFzZSksCiAgICAjbGFtYmRhLjFzZS5tZWQ9bWVkaWFuKGxhbWJkYS4xc2UpLCAKICAgICNsYW1iZGEuMXNlLmhpZ2g9bGFtYmRhLjFzZStsYW1iZGEuMXNlLnNlbSwKICAgICNsYW1iZGEuMXNlLmxvdz1sYW1iZGEuMXNlLWxhbWJkYS4xc2Uuc2VtLAogICAgbnplcm89bnplcm9bMV0sCiAgICBsYW1iZGE9bGFtYmRhWzFdCiAgKQoKbWV0X3Blcl9tZ19zdW1tYXJ5X2xhbWJkYQpgYGAKCgpwbG90IGl0CmBgYHtyfQptZXRfcGVyX21nX3N1bW1hcnlfY3ZtICU+JQogICNmaWx0ZXIoYWxwaGEhPTApICU+JSAjIHdvcnNlIHRoYW4gZXZlcnl0aGluZyBlbHNlIGFuZCB0aHJvd2luZyB0aGUgcGxvdHMgb2ZmCiAgZ2dwbG90KGFlcyh4PWxvZyhsYW1iZGEpLCB5PSBtZWFuY3ZtLCAgeW1pbj1sb3csIHltYXg9aGlnaCkpICsKICBnZW9tX3JpYmJvbihhbHBoYT0uMjUpICsKICBnZW9tX2xpbmUoYWVzKGNvbG9yPWFzLmNoYXJhY3RlcihhbHBoYSkpKSArCiAgZmFjZXRfd3JhcCh+IGFzLmNoYXJhY3RlcihhbHBoYSkpICsKICAgY29vcmRfY2FydGVzaWFuKHhsaW09KGMoLTUsMCkpKSArCiAgZ2VvbV92bGluZShhZXMoeGludGVyY2VwdD1sb2cobGFtYmRhLm1pbi5tZWFuKSksIGFscGhhPS41LCBkYXRhPW1ldF9wZXJfbWdfc3VtbWFyeV9sYW1iZGEpICsKICBnZW9tX3ZsaW5lKGFlcyh4aW50ZXJjZXB0PWxvZyhsYW1iZGEubWluLmhpZ2gpKSwgYWxwaGE9LjUsIGRhdGE9bWV0X3Blcl9tZ19zdW1tYXJ5X2xhbWJkYSwgY29sb3I9ImJsdWUiKSAKCmBgYAoKCgpTbyBvdmVyYWxsIHRoZXNlIGxvb2sgbW9yZSByZWFzb25hYmxlIHRoYW4gdGhlIExPTyBwbG90LgoKTWFrZSBhIHBsb3Qgb2YgTVNFIGF0IG1pbmltdW0gbGFtYmRhIGZvciBlYWNoIGFscGhhCgpgYGB7cn0KbWV0X3Blcl9tZ19zdW1tYXJ5X2N2bSAlPiUgCiAgZ3JvdXBfYnkoYWxwaGEpICU+JQogIGZpbHRlcihyYW5rKG1lYW5jdm0sIHRpZXMubWV0aG9kID0gImZpcnN0Iik9PTEpICU+JQogIGdncGxvdChhZXMoeD1hbHBoYSx5PW1lYW5jdm0seW1pbj1sb3cseW1heD1oaWdoKSkgKwogIGdlb21fcmliYm9uKGNvbG9yPU5BLCBmaWxsPSJncmF5ODAiKSArCiAgZ2VvbV9saW5lKCkgKwogIGdlb21fcG9pbnQoKQpgYGAKbm90IGEgcGFydGljdWxhciBsYXJnZSBkaWZmZXJlbmNlIHRoZXJlLCBhc2lkZSBmcm9tIDAuMSBhbmQgZXZlbiB0aGVuLCBub3QgdG9vIG11Y2ggYmV0dGVyLgoKUGxvdCB0aGUgbnVtYmVyIG9mIG56ZXJvIGNvZWZmaWNpZW50cwoKYGBge3J9Cm1ldF9wZXJfbWdfc3VtbWFyeV9sYW1iZGEgJT4lCiAgdW5uZXN0KGMobGFtYmRhLCBuemVybykpICU+JQogIGdyb3VwX2J5KGFscGhhKSAlPiUKICBmaWx0ZXIoYWJzKGxhbWJkYS5taW4ubWVhbi1sYW1iZGEpPT1taW4oYWJzKGxhbWJkYS5taW4ubWVhbi1sYW1iZGEpKSAgKSAlPiUKICB1bmdyb3VwKCkgJT4lCgpnZ3Bsb3QoYWVzKHg9YXMuY2hhcmFjdGVyKGFscGhhKSwgeT1uemVybykpICsKICBnZW9tX3BvaW50KCkgKwogIGdndGl0bGUoIk51bWJlciBvZiBub24temVybyBjb2VmZmljZW50cyBhdCBtaW5pbXVtIGxhbWJkYSIpICsKICB5bGltKDAsMjQpCmBgYApPSyBsZXQncyBkbyByZXBlYXRlZCB0ZXN0IHRyYWluIHN0YXJ0aW5nIGZyb20gdGhlc2UgQ1YgbGFtYmRhcwoKYGBge3IsIHdhcm5pbmc9RkFMU0V9Cm11bHRpX3R0IDwtIGZ1bmN0aW9uKGxhbWJkYSwgYWxwaGEsIG49MTAwMDAsIHNhbXBsZV9zaXplPTI0LCB0cmFpbl9zaXplPTIwLCB4LCB5PWxlYWZsZW5ndGgkbGVhZl9hdmdfc3RkKSB7CiAgcHJpbnQobGFtYmRhKQogIHByaW50KGFscGhhKQp0dCA8LQogIHRpYmJsZShydW49MTpuKSAlPiUKICBtdXRhdGUodHJhaW49bWFwKHJ1biwgfiBzYW1wbGUoMTpzYW1wbGVfc2l6ZSwgdHJhaW5fc2l6ZSkpKSAlPiUKICBtdXRhdGUoZml0PW1hcCh0cmFpbiwgfiBnbG1uZXQoeD14Wy4sXSwgeT15Wy5dLCBsYW1iZGEgPSBsYW1iZGEsIGFscGhhID0gYWxwaGEgKSkpICU+JQogIAogIG11dGF0ZShwcmVkPW1hcDIoZml0LCB0cmFpbiwgfiBwcmVkaWN0KC54LCBuZXd4ID0geFstLnksXSkpKSAlPiUKICBtdXRhdGUoY29yPW1hcDJfZGJsKHByZWQsIHRyYWluLCB+IGNvcigueCwgeVstLnldKSAgKSkgJT4lCiAgbXV0YXRlKE1TRT1tYXAyX2RibChwcmVkLCB0cmFpbiwgfiBtZWFuKCh5Wy0ueV0gLSAueCleMikpKSAlPiUKICBzdW1tYXJpemUoCiAgICBudW1fbmE9c3VtKGlzLm5hKGNvcikpLCAKICAgIG51bV9sdF8wPXN1bShjb3I8PTAsIG5hLnJtPVRSVUUpLAogICAgYXZnX2Nvcj1tZWFuKGNvciwgbmEucm09VFJVRSksCiAgICBhdmdfTVNFPW1lYW4oTVNFKSkKdHQKfQoKcGVyX21nX2ZpdF90ZXN0X3RyYWluIDwtIG1ldF9wZXJfbWdfc3VtbWFyeV9sYW1iZGEgJT4lIAogIHNlbGVjdChhbHBoYSwgbGFtYmRhLm1pbi5tZWFuKQoKcGVyX21nX2ZpdF90ZXN0X3RyYWluIDwtIG1ldF9wZXJfbWdfbXVsdGlDViAlPiUKICBmaWx0ZXIocnVuPT0xKSAlPiUKICBzZWxlY3QoYWxwaGEsIGZpdCkgJT4lCiAgcmlnaHRfam9pbihwZXJfbWdfZml0X3Rlc3RfdHJhaW4pCgpwZXJfbWdfZml0X3Rlc3RfdHJhaW4gPC0gcGVyX21nX2ZpdF90ZXN0X3RyYWluICU+JQogIG11dGF0ZShwcmVkX2Z1bGw9bWFwMihmaXQsIGxhbWJkYS5taW4ubWVhbiwgfiBwcmVkaWN0KC54LCBzPS55LCBuZXd4PW1ldF9wZXJfbWcuUENzKSksCiAgICAgICAgIGZ1bGxfUj1tYXBfZGJsKHByZWRfZnVsbCwgfiBjb3IoLngsIGxlYWZsZW5ndGgkbGVhZl9hdmdfc3RkKSksCiAgICAgICAgIGZ1bGxfTVNFPW1hcF9kYmwocHJlZF9mdWxsLCB+IG1lYW4oKGxlYWZsZW5ndGgkbGVhZl9hdmdfc3RkLS54KV4yKSkpICU+JQogIAogIG11dGF0ZSh0dD1tYXAyKGxhbWJkYS5taW4ubWVhbiwgYWxwaGEsIH4gbXVsdGlfdHQobGFtYmRhPS54LCBhbHBoYT0ueSwgeD1tZXRfcGVyX21nLlBDcykpKQoKCgoocGVyX21nX2ZpdF90ZXN0X3RyYWluIDwtIHBlcl9tZ19maXRfdGVzdF90cmFpbiAlPiUgdW5uZXN0KHR0KSkKYGBgCgpgYGB7cn0KcGVyX21nX2ZpdF90ZXN0X3RyYWluICU+JQogIGdncGxvdChhZXMoeD1hbHBoYSkpICsKICBnZW9tX2xpbmUoYWVzKHk9YXZnX2NvciksIGNvbG9yPSJyZWQiKSArCiAgZ2VvbV9wb2ludChhZXMoeT1hdmdfY29yKSwgY29sb3I9InJlZCIpICsKICBnZW9tX2xpbmUoYWVzKHk9YXZnX01TRSksIGNvbG9yPSJibHVlIikgKwogIGdlb21fcG9pbnQoYWVzKHk9YXZnX01TRSksIGNvbG9yPSJibHVlIikKYGBgCmFsbCBhcmUgc2ltaWxhciBmcm9tIDAuMiBvbndhcmRzCgojIyBsb29rIGF0IGZpdDoKCmBgYHtyfQphbHBoYV9wZXJfbWcgPC0gLjIKCmJlc3RfcGVyX21nIDwtIHBlcl9tZ19maXRfdGVzdF90cmFpbiAlPiUgZmlsdGVyKGFscGhhID09IGFscGhhX3Blcl9tZykgCmJlc3RfcGVyX21nX2ZpdCA8LSBiZXN0X3Blcl9tZyRmaXRbWzFdXQpiZXN0X3Blcl9tZ19sYW1iZGEgPC0gYmVzdF9wZXJfbWckbGFtYmRhLm1pbi5tZWFuCgpwZXJfbWdfY29lZi50YiA8LSBjb2VmKGJlc3RfcGVyX21nX2ZpdCwgcz1iZXN0X3Blcl9tZ19sYW1iZGEpICU+JSAKICBhcy5tYXRyaXgoKSAlPiUgYXMuZGF0YS5mcmFtZSgpICU+JSAKICByb3duYW1lc190b19jb2x1bW4odmFyPSJQQyIpICU+JQogIHJlbmFtZShiZXRhPWAxYCkKICAKcGVyX21nX2NvZWYudGIgJT4lIGZpbHRlcihiZXRhIT0wKSAlPiUgYXJyYW5nZShiZXRhKQoKYGBgCgpwcmVkIGFuZCBvYnMKYGBge3J9CnBsb3QobGVhZmxlbmd0aCRsZWFmX2F2Z19zdGQsIGJlc3RfcGVyX21nJHByZWRfZnVsbFtbMV1dKQpjb3IudGVzdChsZWFmbGVuZ3RoJGxlYWZfYXZnX3N0ZCwgYmVzdF9wZXJfbWckcHJlZF9mdWxsW1sxXV0pICMuNjUKYmVzdF9wZXJfbWckZnVsbF9NU0UKYGBgCgojIyBQZXJjZW50IHZhcmlhbmNlIGV4cGxhaW5lZAoKYGBge3J9CnBlcl9tZ192YXJzIDwtIHBlcl9tZ19jb2VmLnRiICU+JSAKICBmaWx0ZXIoYmV0YSAhPTAsIFBDIT0iKEludGVyY2VwdCkiKSAlPiUKICBwdWxsKFBDKSAlPiUgYygibGVhZl9hdmdfc3RkIiwgLikKCnBlcl9tZ19yZWxpbXAgPC0gbGVhZmxlbmd0aCAlPiUgc2VsZWN0KGxlYWZfYXZnX3N0ZCkgJT4lIGNiaW5kKG1ldF9wZXJfbWcuUENzKSAlPiUgYXMuZGF0YS5mcmFtZSgpICU+JSBkcGx5cjo6c2VsZWN0KGFsbF9vZihwZXJfbWdfdmFycykpICU+JQogIGNhbGMucmVsaW1wKCkgCgpwZXJfbWdfY29lZi50YiA8LSBwZXJfbWdfcmVsaW1wQGxtZyAlPiUgYXMubWF0cml4KCkgJT4lIGFzLmRhdGEuZnJhbWUoKSAlPiUKICByb3duYW1lc190b19jb2x1bW4oIlBDIikgJT4lCiAgcmVuYW1lKFByb3BWYXJfbWV0X3Blcl9tZz1WMSkgJT4lCiAgZnVsbF9qb2luKHBlcl9tZ19jb2VmLnRiKSAlPiUKICBhcnJhbmdlKGRlc2MoUHJvcFZhcl9tZXRfcGVyX21nKSkKCnBlcl9tZ19jb2VmLnRiCgpgYGAKCkNoZWNrb3V0IHRoZSByb3RhdGlvbnMuICAKCmBgYHtyfQptZXRfcGVyX21nX3JvdGF0aW9uX291dCA8LSBtZXRfcGVyX21nLlBDX3JvdGF0aW9uICU+JSAKICBwaXZvdF9sb25nZXIoLW1ldGFib2xpdGUsIG5hbWVzX3RvPSJQQyIsIHZhbHVlc190bz0ibG9hZGluZyIpICU+JQogIGZpbHRlcihQQyAlaW4lIGZpbHRlcihwZXJfbWdfY29lZi50YiwgYmV0YSE9MCkkUEMgKSAlPiUKICBncm91cF9ieShQQykgJT4lCiAgICBmaWx0ZXIoIXN0cl9kZXRlY3QobWV0YWJvbGl0ZSwiLioobGVhZnxyb290KV9bMC05XSokIikpICU+JQogIGZpbHRlcihhYnMobG9hZGluZykgPj0gMC4wNSkgJT4lCiAgbGVmdF9qb2luKHBlcl9tZ19jb2VmLnRiLCBieT0iUEMiKSAlPiUKICBhcnJhbmdlKGRlc2MoYWJzKGJldGEpKSwgZGVzYyhhYnMobG9hZGluZykpKSAlPiUKICBtdXRhdGUob3JnYW49aWZlbHNlKHN0cl9kZXRlY3QobWV0YWJvbGl0ZSwgIl9sZWFmXyIpLCAibGVhZiIsICJyb290IiksCiAgICAgICAgIHRyYW5zZm9ybWF0aW9uPSJub3JtYWxpemVkIiwKICAgICAgICAgbWV0YWJvbGl0ZT1zdHJfcmVtb3ZlKG1ldGFib2xpdGUsICJtZXRfcGVyX21nXyhyb290fGxlYWYpXyIpLAogICAgICAgICBtZXRhYm9saXRlX2VmZmVjdF9vbl9sZWFmPWlmZWxzZShiZXRhKmxvYWRpbmc+MCwgImluY3JlYXNlIiwgImRlY3JlYXNlIikpCm1ldF9wZXJfbWdfcm90YXRpb25fb3V0ICU+JSAgd3JpdGVfY3N2KCIuLi9vdXRwdXQvTGVhZl9hc3NvY2lhdGVkX21ldGFib2xpdGVzX25vcm1hbGl6ZWRfTk9CTEFOSy5jc3YiKQptZXRfcGVyX21nX3JvdGF0aW9uX291dApgYGAKCiMgbm9uLW5vcm1hemxpemVkCgojIyBtdWx0aSBDVgoKRml0IDEwMSBDVnMgZm9yIGVhY2ggb2YgMTEgYWxwaGFzCmBgYHtyfQpzZXQuc2VlZCgxMjQ1KQoKZm9sZHMgPC0gdGliYmxlKHJ1bj0xOjEwMSkgJT4lIAogIG11dGF0ZShmb2xkcz1tYXAocnVuLCB+IHNhbXBsZShyZXAoMTo2LDQpKSkpCgpzeXN0ZW0udGltZSAobWV0X2FtdF9tdWx0aUNWIDwtIGV4cGFuZF9ncmlkKHJ1bj0xOjEwMCwgYWxwaGE9cm91bmQoc2VxKDAsMSwuMSksMSkpICU+JQogICAgICAgICAgICAgICBsZWZ0X2pvaW4oZm9sZHMsIGJ5PSJydW4iKSAlPiUKICAgICAgICAgICAgICAgbXV0YXRlKGZpdD1tYXAyKGZvbGRzLCBhbHBoYSwgfiBjdi5nbG1uZXQoeD1tZXRfYW10LlBDcywgeT1sZWFmbGVuZ3RoJGxlYWZfYXZnX3N0ZCwgZm9sZGlkID0gLngsIGFscGhhPS55CiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICkpKQogICAgICAgICAgICAgIywgbGFtYmRhPWV4cChzZXEoLTUsMCxsZW5ndGgub3V0ID0gNTApKSApKSkKKSAjMTAwIHNlY29uZHMKCmhlYWQobWV0X2FtdF9tdWx0aUNWKQpgYGAKCmZvciBlYWNoIGZpdCwgcHVsbCBvdXQgdGhlIG1lYW4gY3YgZXJyb3IsIGxhbWJkYSwgbWluIGxhbWJkYSwgYW5kIDFzZSBsYW1iZGEgCmBgYHtyfQptZXRfYW10X211bHRpQ1YgPC0gbWV0X2FtdF9tdWx0aUNWICU+JQogIG11dGF0ZShjdm09bWFwKGZpdCwgbWFncml0dHI6OmV4dHJhY3QoImN2bSIpKSwKICAgICAgICAgbGFtYmRhPW1hcChmaXQsIG1hZ3JpdHRyOjpleHRyYWN0KCJsYW1iZGEiKSksCiAgICAgICAgIGxhbWJkYS5taW49bWFwX2RibChmaXQsIG1hZ3JpdHRyOjpleHRyYWN0KCJsYW1iZGEubWluIiApKSwKICAgICAgICAgbGFtYmRhLjFzZT1tYXBfZGJsKGZpdCwgbWFncml0dHI6OmV4dHJhY3QoImxhbWJkYS4xc2UiKSksCiAgICAgICAgIG56ZXJvPW1hcChmaXQsIG1hZ3JpdHRyOjpleHRyYWN0KCJuemVybyIpKQogICkKCmhlYWQobWV0X2FtdF9tdWx0aUNWKQpgYGAKCgpub3cgY2FsY3VsYXRlIHRoZSBtZWFuIGFuZCBzZW0gb2YgY3ZtIGFuZCBtaW4sMXNlIGxhYm1kYXMuICBUaGVzZSBuZWVkIHRvIGJlIGRvbmUgc2VwYXJhdGVseSBiZWNhdXNlIG9mIHRoZSB3YXkgdGhlIGdyb3VwaW5nIHdvcmtzCmBgYHtyfQptZXRfYW10X3N1bW1hcnlfY3ZtIDwtIG1ldF9hbXRfbXVsdGlDViAlPiUgZHBseXI6OnNlbGVjdCgtZml0LCAtZm9sZHMpICU+JSAKICB1bm5lc3QoYyhjdm0sIGxhbWJkYSkpICU+JQogIGdyb3VwX2J5KGFscGhhLCBsYW1iZGEpICU+JQogIHN1bW1hcml6ZShtZWFuY3ZtPW1lYW4oY3ZtKSwgc2VtPXNkKGN2bSkvc3FydChuKCkpLCBoaWdoPW1lYW5jdm0rc2VtLCBsb3c9bWVhbmN2bS1zZW0pCgptZXRfYW10X3N1bW1hcnlfY3ZtCmBgYAoKYGBge3J9Cm1ldF9hbXRfc3VtbWFyeV9sYW1iZGEgPC0gbWV0X2FtdF9tdWx0aUNWICU+JSBkcGx5cjo6c2VsZWN0KC1maXQsIC1mb2xkcywgLWN2bSkgJT4lIAogIGdyb3VwX2J5KGFscGhhKSAlPiUKICBzdW1tYXJpemUoCiAgICBsYW1iZGEubWluLnNkPXNkKGxhbWJkYS5taW4pLCAKICAgIGxhbWJkYS5taW4ubWVhbj1tZWFuKGxhbWJkYS5taW4pLAogICAgI2xhbWJkYS5taW4ubWVkPW1lZGlhbihsYW1iZGEubWluKSwgCiAgICBsYW1iZGEubWluLmhpZ2g9bGFtYmRhLm1pbi5tZWFuK2xhbWJkYS5taW4uc2QsCiAgICAjbGFtYmRhLm1pbi5sb3c9bGFtYmRhLm1pbi5tZWFuLWxhbWJkYS5taW4uc2VtLAogICAgI2xhbWJkYS4xc2Uuc2VtPXNkKGxhbWJkYS4xc2UpL3NxcnQobigpKSwgCiAgICBsYW1iZGEuMXNlLm1lYW49bWVhbihsYW1iZGEuMXNlKSwKICAgICNsYW1iZGEuMXNlLm1lZD1tZWRpYW4obGFtYmRhLjFzZSksIAogICAgI2xhbWJkYS4xc2UuaGlnaD1sYW1iZGEuMXNlK2xhbWJkYS4xc2Uuc2VtLAogICAgI2xhbWJkYS4xc2UubG93PWxhbWJkYS4xc2UtbGFtYmRhLjFzZS5zZW0sCiAgICBuemVybz1uemVyb1sxXSwKICAgIGxhbWJkYT1sYW1iZGFbMV0KICApCgptZXRfYW10X3N1bW1hcnlfbGFtYmRhCmBgYAoKCnBsb3QgaXQKYGBge3J9Cm1ldF9hbXRfc3VtbWFyeV9jdm0gJT4lCiAgI2ZpbHRlcihhbHBoYSE9MCkgJT4lICMgd29yc2UgdGhhbiBldmVyeXRoaW5nIGVsc2UgYW5kIHRocm93aW5nIHRoZSBwbG90cyBvZmYKICBnZ3Bsb3QoYWVzKHg9bG9nKGxhbWJkYSksIHk9IG1lYW5jdm0sICB5bWluPWxvdywgeW1heD1oaWdoKSkgKwogIGdlb21fcmliYm9uKGFscGhhPS4yNSkgKwogIGdlb21fbGluZShhZXMoY29sb3I9YXMuY2hhcmFjdGVyKGFscGhhKSkpICsKICBmYWNldF93cmFwKH4gYXMuY2hhcmFjdGVyKGFscGhhKSkgKwogICBjb29yZF9jYXJ0ZXNpYW4oeGxpbT0oYygtNSwwKSkpICsKICBnZW9tX3ZsaW5lKGFlcyh4aW50ZXJjZXB0PWxvZyhsYW1iZGEubWluLm1lYW4pKSwgYWxwaGE9LjUsIGRhdGE9bWV0X2FtdF9zdW1tYXJ5X2xhbWJkYSkgKwogIGdlb21fdmxpbmUoYWVzKHhpbnRlcmNlcHQ9bG9nKGxhbWJkYS5taW4uaGlnaCkpLCBhbHBoYT0uNSwgZGF0YT1tZXRfYW10X3N1bW1hcnlfbGFtYmRhLCBjb2xvcj0iYmx1ZSIpIAoKYGBgCgoKTWFrZSBhIHBsb3Qgb2YgTVNFIGF0IG1pbmltdW0gbGFtYmRhIGZvciBlYWNoIGFscGhhCgpgYGB7cn0KbWV0X2FtdF9zdW1tYXJ5X2N2bSAlPiUgCiAgZ3JvdXBfYnkoYWxwaGEpICU+JQogIGZpbHRlcihyYW5rKG1lYW5jdm0sIHRpZXMubWV0aG9kID0gImZpcnN0Iik9PTEpICU+JQogIGdncGxvdChhZXMoeD1hbHBoYSx5PW1lYW5jdm0seW1pbj1sb3cseW1heD1oaWdoKSkgKwogIGdlb21fcmliYm9uKGNvbG9yPU5BLCBmaWxsPSJncmF5ODAiKSArCiAgZ2VvbV9saW5lKCkgKwogIGdlb21fcG9pbnQoKQpgYGAKZGVjcmVhc2luZyB0aHJvdWdob3V0LCBhbHRob3VnaCBkaWZmZXJlbmNlcyBhcmUgcHJldHR5IHNtYWxsCgpQbG90IHRoZSBudW1iZXIgb2Ygbnplcm8gY29lZmZpY2llbnRzCgpgYGB7cn0KbWV0X2FtdF9zdW1tYXJ5X2xhbWJkYSAlPiUKICB1bm5lc3QoYyhsYW1iZGEsIG56ZXJvKSkgJT4lCiAgZ3JvdXBfYnkoYWxwaGEpICU+JQogIGZpbHRlcihhYnMobGFtYmRhLm1pbi5tZWFuLWxhbWJkYSk9PW1pbihhYnMobGFtYmRhLm1pbi5tZWFuLWxhbWJkYSkpICApICU+JQogIHVuZ3JvdXAoKSAlPiUKCmdncGxvdChhZXMoeD1hcy5jaGFyYWN0ZXIoYWxwaGEpLCB5PW56ZXJvKSkgKwogIGdlb21fcG9pbnQoKSArCiAgZ2d0aXRsZSgiTnVtYmVyIG9mIG5vbi16ZXJvIGNvZWZmaWNlbnRzIGF0IG1pbmltdW0gbGFtYmRhIikgKwogIHlsaW0oMCwzNikKYGBgCk9LIGxldCdzIGRvIHJlcGVhdGVkIHRlc3QgdHJhaW4gc3RhcnRpbmcgZnJvbSB0aGVzZSBDViBsYW1iZGFzCgpgYGB7cn0KbXVsdGlfdHQgPC0gZnVuY3Rpb24obGFtYmRhLCBhbHBoYSwgbj0xMDAwMCwgc2FtcGxlX3NpemU9MjQsIHRyYWluX3NpemU9MjAsIHgsIHk9bGVhZmxlbmd0aCRsZWFmX2F2Z19zdGQpIHsKICBwcmludChsYW1iZGEpCiAgcHJpbnQoYWxwaGEpCnR0IDwtCiAgdGliYmxlKHJ1bj0xOm4pICU+JQogIG11dGF0ZSh0cmFpbj1tYXAocnVuLCB+IHNhbXBsZSgxOnNhbXBsZV9zaXplLCB0cmFpbl9zaXplKSkpICU+JQogIG11dGF0ZShmaXQ9bWFwKHRyYWluLCB+IGdsbW5ldCh4PXhbLixdLCB5PXlbLl0sIGxhbWJkYSA9IGxhbWJkYSwgYWxwaGEgPSBhbHBoYSApKSkgJT4lCiAgCiAgbXV0YXRlKHByZWQ9bWFwMihmaXQsIHRyYWluLCB+IHByZWRpY3QoLngsIG5ld3ggPSB4Wy0ueSxdKSkpICU+JQogIG11dGF0ZShjb3I9bWFwMl9kYmwocHJlZCwgdHJhaW4sIH4gY29yKC54LCB5Wy0ueV0pICApKSAlPiUKICBtdXRhdGUoTVNFPW1hcDJfZGJsKHByZWQsIHRyYWluLCB+IG1lYW4oKHlbLS55XSAtIC54KV4yKSkpICU+JQogIHN1bW1hcml6ZSgKICAgIG51bV9uYT1zdW0oaXMubmEoY29yKSksIAogICAgbnVtX2x0XzA9c3VtKGNvcjw9MCwgbmEucm09VFJVRSksCiAgICBhdmdfY29yPW1lYW4oY29yLCBuYS5ybT1UUlVFKSwKICAgIGF2Z19NU0U9bWVhbihNU0UpKQp0dAp9CgphbXRfZml0X3Rlc3RfdHJhaW4gPC0gbWV0X2FtdF9zdW1tYXJ5X2xhbWJkYSAlPiUgCiAgc2VsZWN0KGFscGhhLCBsYW1iZGEubWluLm1lYW4pCgphbXRfZml0X3Rlc3RfdHJhaW4gPC0gbWV0X2FtdF9tdWx0aUNWICU+JQogIGZpbHRlcihydW49PTEpICU+JQogIHNlbGVjdChhbHBoYSwgZml0KSAlPiUKICByaWdodF9qb2luKGFtdF9maXRfdGVzdF90cmFpbikKCmFtdF9maXRfdGVzdF90cmFpbiA8LSBhbXRfZml0X3Rlc3RfdHJhaW4gJT4lCiAgbXV0YXRlKHByZWRfZnVsbD1tYXAyKGZpdCwgbGFtYmRhLm1pbi5tZWFuLCB+IHByZWRpY3QoLngsIHM9LnksIG5ld3g9bWV0X2FtdC5QQ3MpKSwKICAgICAgICAgZnVsbF9SPW1hcF9kYmwocHJlZF9mdWxsLCB+IGNvcigueCwgbGVhZmxlbmd0aCRsZWFmX2F2Z19zdGQpKSwKICAgICAgICAgZnVsbF9NU0U9bWFwX2RibChwcmVkX2Z1bGwsIH4gbWVhbigobGVhZmxlbmd0aCRsZWFmX2F2Z19zdGQtLngpXjIpKSkgJT4lCiAgCiAgbXV0YXRlKHR0PW1hcDIobGFtYmRhLm1pbi5tZWFuLCBhbHBoYSwgfiBtdWx0aV90dChsYW1iZGE9LngsIGFscGhhPS55LCB4PW1ldF9hbXQuUENzKSkpCgoKCihhbXRfZml0X3Rlc3RfdHJhaW4gPC0gYW10X2ZpdF90ZXN0X3RyYWluICU+JSB1bm5lc3QodHQpKQpgYGAKCmBgYHtyfQphbXRfZml0X3Rlc3RfdHJhaW4gJT4lCiAgZ2dwbG90KGFlcyh4PWFscGhhKSkgKwogIGdlb21fbGluZShhZXMoeT1hdmdfY29yKSwgY29sb3I9InJlZCIpICsKICBnZW9tX3BvaW50KGFlcyh5PWF2Z19jb3IpLCBjb2xvcj0icmVkIikgKwogIGdlb21fbGluZShhZXMoeT1hdmdfTVNFKSwgY29sb3I9ImJsdWUiKSArCiAgZ2VvbV9wb2ludChhZXMoeT1hdmdfTVNFKSwgY29sb3I9ImJsdWUiKQpgYGAKbm90IGEgYmlnIGRpZmZlcmVuY2UgYWZ0ZXIgMC4yCgojIyBsb29rIGF0IGZpdDoKCmBgYHtyfQphbHBoYV9hbXQgPC0gLjcKCmJlc3RfYW10IDwtIGFtdF9maXRfdGVzdF90cmFpbiAlPiUgZmlsdGVyKGFscGhhID09IGFscGhhX2FtdCkgCmJlc3RfYW10X2ZpdCA8LSBiZXN0X2FtdCRmaXRbWzFdXQpiZXN0X2FtdF9sYW1iZGEgPC0gYmVzdF9hbXQkbGFtYmRhLm1pbi5tZWFuCgphbXRfY29lZi50YiA8LSBjb2VmKGJlc3RfYW10X2ZpdCwgcz1iZXN0X2FtdF9sYW1iZGEpICU+JSAKICBhcy5tYXRyaXgoKSAlPiUgYXMuZGF0YS5mcmFtZSgpICU+JSAKICByb3duYW1lc190b19jb2x1bW4odmFyPSJQQyIpICU+JQogIHJlbmFtZShiZXRhPWAxYCkKICAKYW10X2NvZWYudGIgJT4lIGZpbHRlcihiZXRhIT0wKSAlPiUgYXJyYW5nZShiZXRhKQoKYGBgCgpwcmVkIGFuZCBvYnMKYGBge3J9CnBsb3QobGVhZmxlbmd0aCRsZWFmX2F2Z19zdGQsIGJlc3RfYW10JHByZWRfZnVsbFtbMV1dKQpjb3IudGVzdChsZWFmbGVuZ3RoJGxlYWZfYXZnX3N0ZCwgYmVzdF9hbXQkcHJlZF9mdWxsW1sxXV0pICMuNjQKYmVzdF9hbXQkZnVsbF9NU0UKYGBgCgojIyBQZXJjZW50IHZhcmlhbmNlIGV4cGxhaW5lZAoKYGBge3J9CmFtdF92YXJzIDwtIGFtdF9jb2VmLnRiICU+JSAKICBmaWx0ZXIoYmV0YSAhPTAsIFBDIT0iKEludGVyY2VwdCkiKSAlPiUKICBwdWxsKFBDKSAlPiUgYygibGVhZl9hdmdfc3RkIiwgLikKCiMgYmVjYXVzZSB0aGVyZSBpcyBvbmx5IG9uZSBleHBsYW5hdG9yeSB2YXJpYWJsZSwgd2UgZG9uJ3QgbmVlZCB0byAoYW5kIGNhbid0KSB1c2UgcmVsaW1wCmFtdF9yZWxpbXAgPC0gbGVhZmxlbmd0aCAlPiUgc2VsZWN0KGxlYWZfYXZnX3N0ZCkgJT4lIGNiaW5kKG1ldF9hbXQuUENzKSAlPiUgYXMuZGF0YS5mcmFtZSgpICU+JSBkcGx5cjo6c2VsZWN0KGFsbF9vZihhbXRfdmFycykpICU+JSBjb3IoKSAlPiUgbWFncml0dHI6OmV4dHJhY3QoMSwyKSAlPiUgbWFncml0dHI6Om11bHRpcGx5X2J5KC4sLikKYW10X3JlbGltcCA8LSB0aWJibGUoUEM9c3RyX3N1YnNldChhbXRfdmFycywgIlBDIiksCiAgICAgICAgICAgICAgICAgICAgIFByb3BWYXJfbWV0X2FtdD1hbXRfcmVsaW1wKQoKCmFtdF9jb2VmLnRiIDwtIGFtdF9yZWxpbXAgJT4lIAogIGZ1bGxfam9pbihhbXRfY29lZi50YikgJT4lCiAgYXJyYW5nZShkZXNjKFByb3BWYXJfbWV0X2FtdCkpCgphbXRfY29lZi50YgoKYGBgCgpDaGVja291dCB0aGUgcm90YXRpb25zLiAgCgpgYGB7cn0KbWV0X2FtdF9yb3RhdGlvbl9vdXQgPC0gbWV0X2FtdC5QQ19yb3RhdGlvbiAlPiUgCiAgcGl2b3RfbG9uZ2VyKC1tZXRhYm9saXRlLCBuYW1lc190bz0iUEMiLCB2YWx1ZXNfdG89ImxvYWRpbmciKSAlPiUKICBmaWx0ZXIoUEMgJWluJSBmaWx0ZXIoYW10X2NvZWYudGIsIGJldGEhPTApJFBDICkgJT4lCiAgZ3JvdXBfYnkoUEMpICU+JQogICAgZmlsdGVyKCFzdHJfZGV0ZWN0KG1ldGFib2xpdGUsIi4qKGxlYWZ8cm9vdClfWzAtOV0qJCIpKSAlPiUKICBmaWx0ZXIoYWJzKGxvYWRpbmcpID49IDAuMDUpICU+JQogIGxlZnRfam9pbihhbXRfY29lZi50YiwgYnk9IlBDIikgJT4lCiAgYXJyYW5nZShkZXNjKGFicyhiZXRhKSksIGRlc2MoYWJzKGxvYWRpbmcpKSkgJT4lCiAgbXV0YXRlKG9yZ2FuPWlmZWxzZShzdHJfZGV0ZWN0KG1ldGFib2xpdGUsICJfbGVhZl8iKSwgImxlYWYiLCAicm9vdCIpLAogICAgICAgICB0cmFuc2Zvcm1hdGlvbj0icmF3IiwKICAgICAgICAgbWV0YWJvbGl0ZT1zdHJfcmVtb3ZlKG1ldGFib2xpdGUsICJtZXRfYW10Xyhyb290fGxlYWYpXyIpLAogICAgICAgICBtZXRhYm9saXRlX2VmZmVjdF9vbl9sZWFmPWlmZWxzZShiZXRhKmxvYWRpbmc+MCwgImluY3JlYXNlIiwgImRlY3JlYXNlIikpCm1ldF9hbXRfcm90YXRpb25fb3V0ICU+JSAgd3JpdGVfY3N2KCIuLi9vdXRwdXQvTGVhZl9hc3NvY2lhdGVkX21ldGFib2xpdGVzX3Jhd19OT0JMQU5LLmNzdiIpCgptZXRfYW10X3JvdGF0aW9uX291dApgYGA=